<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/>
<meta http-equiv="X-UA-Compatible" content="IE=9"/>
<title>AMBHAS: /home/tomer/svn/ambhas/ambhas/rain_disagg.py Source File</title>

<link href="tabs.css" rel="stylesheet" type="text/css"/>
<link href="doxygen.css" rel="stylesheet" type="text/css" />
<link href="navtree.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="jquery.js"></script>
<script type="text/javascript" src="resize.js"></script>
<script type="text/javascript" src="navtree.js"></script>
<script type="text/javascript">
  $(document).ready(initResizable);
</script>
<link href="search/search.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="search/search.js"></script>
<script type="text/javascript">
  $(document).ready(function() { searchBox.OnSelectItem(0); });
</script>

</head>
<body>
<div id="top"><!-- do not remove this div! -->


<div id="titlearea">
<table cellspacing="0" cellpadding="0">
 <tbody>
 <tr style="height: 56px;">
  
  
  <td style="padding-left: 0.5em;">
   <div id="projectname">AMBHAS
   
   </div>
   
  </td>
  
  
  
   
   <td>        <div id="MSearchBox" class="MSearchBoxInactive">
        <span class="left">
          <img id="MSearchSelect" src="search/mag_sel.png"
               onmouseover="return searchBox.OnSearchSelectShow()"
               onmouseout="return searchBox.OnSearchSelectHide()"
               alt=""/>
          <input type="text" id="MSearchField" value="Search" accesskey="S"
               onfocus="searchBox.OnSearchFieldFocus(true)" 
               onblur="searchBox.OnSearchFieldFocus(false)" 
               onkeyup="searchBox.OnSearchFieldChange(event)"/>
          </span><span class="right">
            <a id="MSearchClose" href="javascript:searchBox.CloseResultsWindow()"><img id="MSearchCloseImg" border="0" src="search/close.png" alt=""/></a>
          </span>
        </div>
</td>
   
  
 </tr>
 </tbody>
</table>
</div>

<!-- Generated by Doxygen 1.7.6.1 -->
<script type="text/javascript">
var searchBox = new SearchBox("searchBox", "search",false,'Search');
</script>
</div>
<div id="side-nav" class="ui-resizable side-nav-resizable">
  <div id="nav-tree">
    <div id="nav-tree-contents">
    </div>
  </div>
  <div id="splitbar" style="-moz-user-select:none;" 
       class="ui-resizable-handle">
  </div>
</div>
<script type="text/javascript">
  initNavTree('rain__disagg_8py.html','');
</script>
<div id="doc-content">
<div class="header">
  <div class="headertitle">
<div class="title">rain_disagg.py</div>  </div>
</div><!--header-->
<div class="contents">
<a href="rain__disagg_8py.html">Go to the documentation of this file.</a><div class="fragment"><pre class="fragment"><a name="l00001"></a><a class="code" href="namespaceambhas_1_1rain__disagg.html">00001</a> <span class="comment"># -*- coding: utf-8 -*-</span>
<a name="l00002"></a>00002 <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00003"></a>00003 <span class="stringliteral">Created on Tue May 24 18:07:28 2011</span>
<a name="l00004"></a>00004 <span class="stringliteral"></span>
<a name="l00005"></a>00005 <span class="stringliteral">@author: Sat Kumar Tomer</span>
<a name="l00006"></a>00006 <span class="stringliteral">@website: www.ambhas.com</span>
<a name="l00007"></a>00007 <span class="stringliteral">@email: satkumartomer@gmail.com</span>
<a name="l00008"></a>00008 <span class="stringliteral">&quot;&quot;&quot;</span>
<a name="l00009"></a>00009 
<a name="l00010"></a>00010 <span class="comment"># import required libraries</span>
<a name="l00011"></a>00011 <span class="keyword">from</span> __future__ <span class="keyword">import</span> division
<a name="l00012"></a>00012 <span class="keyword">import</span> numpy <span class="keyword">as</span> np
<a name="l00013"></a>00013 <span class="keyword">from</span> ambhas.errlib <span class="keyword">import</span> rmse
<a name="l00014"></a>00014 <span class="keyword">from</span> scipy.optimize <span class="keyword">import</span> fmin
<a name="l00015"></a>00015 <span class="keyword">from</span> scipy.stats <span class="keyword">import</span> poisson
<a name="l00016"></a>00016 
<a name="l00017"></a><a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html">00017</a> <span class="keyword">class </span><a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html">RainDisagg</a>:
<a name="l00018"></a>00018     
<a name="l00019"></a>00019     
<a name="l00020"></a><a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html#ad31d5b9ac3e571d6f1e199f6b9a63831">00020</a>     <span class="keyword">def </span><a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html#ab38afc38eb480247f09d0e8fee937876">__init__</a>(self,rf):
<a name="l00021"></a>00021         self.<a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html#a7268936d7403b1b8d40b0f8060c3fad6">rf</a> = rf
<a name="l00022"></a>00022             
<a name="l00023"></a>00023         <span class="comment"># calculate the length of rainfall series</span>
<a name="l00024"></a>00024         len_rf = len(rf)
<a name="l00025"></a>00025                   
<a name="l00026"></a>00026         <span class="comment"># take only that much length of rainfall series which is multiplication of 32 </span>
<a name="l00027"></a>00027         rf_1 = rf[0:len_rf-np.mod(len_rf,32)]
<a name="l00028"></a>00028         
<a name="l00029"></a>00029         <span class="comment"># summing rainfall for 2 days</span>
<a name="l00030"></a>00030         rf_2 = np.sum(rf_1.reshape(-1,2), axis=1)
<a name="l00031"></a>00031         
<a name="l00032"></a>00032         <span class="comment"># summing rainfall for 4 days</span>
<a name="l00033"></a>00033         rf_4 = np.sum(rf_2.reshape(-1,2), axis=1)
<a name="l00034"></a>00034         
<a name="l00035"></a>00035         <span class="comment"># summing rainfall for 8 days</span>
<a name="l00036"></a>00036         rf_8 = np.sum(rf_4.reshape(-1,2), axis=1)
<a name="l00037"></a>00037         
<a name="l00038"></a>00038         <span class="comment"># summing rainfall for 16 days</span>
<a name="l00039"></a>00039         rf_16 = np.sum(rf_8.reshape(-1,2), axis=1)
<a name="l00040"></a>00040         
<a name="l00041"></a>00041         <span class="comment"># summing rainfall for 32 days</span>
<a name="l00042"></a>00042         rf_32 = np.sum(rf_16.reshape(-1,2), axis=1)
<a name="l00043"></a>00043         
<a name="l00044"></a>00044         <span class="comment">#generate moments (q varies from 0 to 5)</span>
<a name="l00045"></a>00045         <span class="comment"># row =&gt; time scale</span>
<a name="l00046"></a>00046         <span class="comment"># column =&gt; moments</span>
<a name="l00047"></a>00047         M1 = np.zeros((6,11))
<a name="l00048"></a>00048         <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(11):
<a name="l00049"></a>00049             M1[0,i] = np.mean(rf_1**(i/2))
<a name="l00050"></a>00050             M1[1,i] = np.mean(rf_2**(i/2))
<a name="l00051"></a>00051             M1[2,i] = np.mean(rf_4**(i/2))
<a name="l00052"></a>00052             M1[3,i] = np.mean(rf_8**(i/2))
<a name="l00053"></a>00053             M1[4,i] = np.mean(rf_16**(i/2))
<a name="l00054"></a>00054             M1[5,i] = np.mean(rf_32**(i/2))
<a name="l00055"></a>00055         
<a name="l00056"></a>00056         self.<a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html#acfce7c0dc9dfadd699f2fe0a10049ca0">M1</a> = M1
<a name="l00057"></a>00057         self.<a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html#a337eed66bd7b4ddb7c4061e5bea495f5">logM</a> = np.log(M1)
<a name="l00058"></a>00058         
<a name="l00059"></a>00059         <span class="comment"># l is inverse of time scale</span>
<a name="l00060"></a>00060         l = [32, 16, 8, 4, 2, 1]
<a name="l00061"></a>00061         self.<a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html#af3b6ec9f7adaa19f9a0dd2871c6876a7">log_lambda</a> = np.log(l)
<a name="l00062"></a>00062         
<a name="l00063"></a>00063         
<a name="l00064"></a>00064         <span class="comment"># calculate the tau from the slope of log(M) vs log(l)</span>
<a name="l00065"></a>00065         tau_obs = np.zeros(10,)
<a name="l00066"></a>00066         <span class="keywordflow">for</span> i <span class="keywordflow">in</span> range(10):
<a name="l00067"></a>00067             tau_obs[i] = -np.polyfit(np.log(l), np.log(M1[:,i+1]),1)[0]
<a name="l00068"></a>00068         self.<a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html#ad31d5b9ac3e571d6f1e199f6b9a63831">tau_obs</a> = tau_obs
<a name="l00069"></a>00069         
<a name="l00070"></a>00070         <span class="comment"># fit the log-poisson distribution</span>
<a name="l00071"></a>00071         self.<a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html#a53add807ea7a117f041c8b7685d9d68b">lp</a> = fmin(self.<a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html#a950a4d1cdae7a50557e568a6084951f1">fun_poisson</a>,np.array([0.4, 0.2]))
<a name="l00072"></a>00072         
<a name="l00073"></a>00073         <span class="comment"># calculate the parameter of log poisson distribution form the parameters</span>
<a name="l00074"></a>00074         <span class="comment"># of tau function</span>
<a name="l00075"></a>00075         self.<a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html#a01807c3ec02381e5a166e0cda1e4a837">A</a> = np.exp(self.<a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html#a53add807ea7a117f041c8b7685d9d68b">lp</a>[0]*(1-self.<a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html#a53add807ea7a117f041c8b7685d9d68b">lp</a>[1]))
<a name="l00076"></a>00076 
<a name="l00077"></a><a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html#aed74cef4625da50546d3d79923120193">00077</a>     <span class="keyword">def </span><a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html#aed74cef4625da50546d3d79923120193">tau_predict</a>(self):
<a name="l00078"></a>00078         q = np.arange(0.5,5.5,0.5)
<a name="l00079"></a>00079         c = abs(self.<a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html#a53add807ea7a117f041c8b7685d9d68b">lp</a>[0])
<a name="l00080"></a>00080         beta = abs(self.<a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html#a53add807ea7a117f041c8b7685d9d68b">lp</a>[1])
<a name="l00081"></a>00081         b = 2
<a name="l00082"></a>00082         tau_pred = q-c*(q*(1-beta)+beta**q-1)/(np.log(b))
<a name="l00083"></a>00083         self.<a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html#ab30ce9b40553b4fcf61f09d3c75eb08b">tau_pred</a> = tau_pred
<a name="l00084"></a>00084         self.<a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html#aad369c07f675afdd3ca6229ff0a12f28">q</a> = q
<a name="l00085"></a>00085 
<a name="l00086"></a>00086     <span class="comment">#define the log-poisson function</span>
<a name="l00087"></a><a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html#a950a4d1cdae7a50557e568a6084951f1">00087</a>     <span class="keyword">def </span><a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html#a950a4d1cdae7a50557e568a6084951f1">fun_poisson</a>(self,par):
<a name="l00088"></a>00088         q = np.arange(0.5,5.5,0.5)
<a name="l00089"></a>00089         c = abs(par[0])
<a name="l00090"></a>00090         beta = abs(par[1])
<a name="l00091"></a>00091         b = 2
<a name="l00092"></a>00092         tau_pred = q-c*(q*(1-beta)+beta**q-1)/(np.log(b))
<a name="l00093"></a>00093         f = <a class="code" href="namespaceambhas_1_1errlib.html#a6644b1f966f2f5f2f2e7be6002ccb225">rmse</a>(tau_pred,self.<a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html#ad31d5b9ac3e571d6f1e199f6b9a63831">tau_obs</a>)
<a name="l00094"></a>00094         <span class="keywordflow">return</span> f
<a name="l00095"></a>00095     
<a name="l00096"></a><a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html#a9b10c013059531b29d1b7e63be79fbb1">00096</a>     <span class="keyword">def </span><a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html#a9b10c013059531b29d1b7e63be79fbb1">disaggregate</a>(self,rf):
<a name="l00097"></a>00097         len_rf = len(rf)        
<a name="l00098"></a>00098         <span class="comment"># generating rainfall from t h to t/2 h</span>
<a name="l00099"></a>00099         rf_pre = np.zeros((1,len_rf*2))
<a name="l00100"></a>00100         <span class="keywordflow">for</span> j <span class="keywordflow">in</span> range(1):
<a name="l00101"></a>00101             <span class="keywordflow">for</span> i <span class="keywordflow">in</span> xrange(0,len_rf*2,2):
<a name="l00102"></a>00102                 W = self.<a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html#a01807c3ec02381e5a166e0cda1e4a837">A</a>*(self.<a class="code" href="classambhas_1_1rain__disagg_1_1RainDisagg.html#a53add807ea7a117f041c8b7685d9d68b">lp</a>[1])**poisson.rvs(1, size=2)
<a name="l00103"></a>00103                 W[W&lt;0] = 1e-6
<a name="l00104"></a>00104                 rf_pre[j,i] = rf[int(i/2)]*W[0]/(W[0]+W[1])
<a name="l00105"></a>00105                 rf_pre[j,i+1] = rf[int(i/2)]*W[1]/(W[0]+W[1])              
<a name="l00106"></a>00106                              
<a name="l00107"></a>00107         
<a name="l00108"></a>00108         rf_pre = np.mean(rf_pre, axis=0)
<a name="l00109"></a>00109         
<a name="l00110"></a>00110         <span class="comment"># rounding up the simulated rainfall to the least count of raingauge </span>
<a name="l00111"></a>00111         <span class="keywordflow">for</span> i <span class="keywordflow">in</span> xrange(0,len_rf*2,2):
<a name="l00112"></a>00112             <span class="keywordflow">if</span> np.mod(rf_pre[i],0.5) !=0:
<a name="l00113"></a>00113                 TB = np.mod(rf_pre[i],0.5)
<a name="l00114"></a>00114             <span class="keywordflow">else</span>:
<a name="l00115"></a>00115                 TB = 0
<a name="l00116"></a>00116             
<a name="l00117"></a>00117             rf_pre[i] -= TB
<a name="l00118"></a>00118             rf_pre[i+1] += TB
<a name="l00119"></a>00119             
<a name="l00120"></a>00120             
<a name="l00121"></a>00121         
<a name="l00122"></a>00122         <span class="keywordflow">return</span> rf_pre
<a name="l00123"></a>00123 
</pre></div></div><!-- contents -->
</div>
<!-- window showing the filter options -->
<div id="MSearchSelectWindow"
     onmouseover="return searchBox.OnSearchSelectShow()"
     onmouseout="return searchBox.OnSearchSelectHide()"
     onkeydown="return searchBox.OnSearchSelectKey(event)">
<a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(0)"><span class="SelectionMark">&#160;</span>All</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(1)"><span class="SelectionMark">&#160;</span>Classes</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(2)"><span class="SelectionMark">&#160;</span>Namespaces</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(3)"><span class="SelectionMark">&#160;</span>Files</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(4)"><span class="SelectionMark">&#160;</span>Functions</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(5)"><span class="SelectionMark">&#160;</span>Variables</a></div>

<!-- iframe showing the search results (closed by default) -->
<div id="MSearchResultsWindow">
<iframe src="javascript:void(0)" frameborder="0" 
        name="MSearchResults" id="MSearchResults">
</iframe>
</div>

  <div id="nav-path" class="navpath">
    <ul>
      <li class="navelem"><a class="el" href="rain__disagg_8py.html">rain_disagg.py</a>      </li>

    <li class="footer">Generated on Sat Jul 21 2012 12:26:08 for AMBHAS by
    <a href="http://www.doxygen.org/index.html">
    <img class="footer" src="doxygen.png" alt="doxygen"/></a> 1.7.6.1 </li>
   </ul>
 </div>


</body>
</html>
